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Abstract 

Numerical challenges inherent in algorithms for computing worst Value-at-Risk in homogeneous 
portfolios are identified and solutions as well as words of warning concerning their implemen¬ 
tation are provided. Furthermore, both conceptual and computational improvements to the 
Rearrangement Algorithm for approximating worst Value-at-Risk for portfolios with arbitrary 
marginal loss distributions are given. In particular, a novel Adaptive Rearrangement Algorithm is 
introduced and investigated. These algorithms are implemented using the R package qrmtools. 
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1 Introduction 

An integral part of Quantitative Risk Management is to analyze the one-period ahead vector of 
losses L = (Li,..., Lrf), where Lj represents the loss (a random variable) associated with a given 
business line or risk type with counterparty j, j € {1,... ,d}, over a fixed time horizon. For financial 
institutions, the aggregated loss 


d 

i=i 

is of particular interest. Under Pillar I of the Basel Accords, financial institutions are required to 
set aside capital to manage market, credit and operational risk. To this end a risk measure p{-) is 
used to map the aggregate position L+ to p{L~^) G IR for obtaining the amount of capital required 
to account for future losses over a predetermined time period. As a risk measure, Value-at-Risk 
(VaRa) has been widely adopted by the financial industry since the mid nineties. It is defined as 
the a-quantile of the distribution function Fy+ of , i.e., 


VaRQ(L^) = F^+(a) = inf{x G IR : F^+(x) > a}, 


12 


for more details. A 


where denotes the quantile function of Fj^+; see Embrechts and Hofert 
well known drawback of VaRQ(L'*') as a risk measure is that VaRa(L'*“) is not necessarily subadditive 
unless L follows an elliptical distribution; see, e.g., Embrechts et al. [^, McNeil et al. [I9| p. 241], 
Embrechts et al. and Hofert and McNeil |16| . 

There are various methods for estimating the marginal loss distributions Fi,... ,Fd of Li,..., L^, 
respectively, but capturing the d-variate dependence structure (i.e., the underlying copula C) 
of L is often more difficult. This is due to the fact that typically not much is known about 
C and estimation often not feasible (e.g., for rare-event losses occurring in different geographic 
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regions). In this work we focus on the case where C is unknown; the case of partial information 
about C, is studied by Bernard et ah and Bernard et ah [^. In our case, one only knows that 
VaRa(L+) G [ VaRa lL"^), VaRQ(L+)] where VaR^ ,fL'*~l and VaRQ(L''') denote the best and the 
worst VaRQ(L+) over all distribution functions of L with marginals Fi,..., Fd, respectively. Note 
that this interval can be wide, but financial firms are interested in computing it (often in high 
dimensions d) to determine their risk capital for L+ within this range. As we show in this work, 
even for small d (and other moderate parameter choices) this can be challenging. 

In particular, we investigate solutions in the homogeneous case (i.e., Fi = ■ ■ ■ = F^) which are 
considered “explicit”; see Embrechts et al. (note that the formulas for both VaR. „,(L'*~) and 
VaRa(E'*') have typographical errors; we correct them below). In the general, inhomogeneous 
case (i.e., not all_K’s necessarily being equal), we consider the Rearrangement Algorithm of 
Embrechts et al. |ll| for computing VaRQ,(L'’“) and VaRQ(L'’'). The presented algorithms with 
conceptual and numerical improvements have been implemented in the R package qrmtools; see 
also the accompanying vignette VaR_bounds which provides further results, numerical investigations, 
diagnostic checks and an application. All the results in this paper can be reproduced with the 
package and vignette (and, obviously, other parameters can be chosen if of interest). For a different 
approach for computing VaR. „,(L~*~) and VaRQ(T“'") not discussed here, see Bernard and McLeish [^. 
In what follows, we focus on the worst VaRQ,(L'*“), i.e., VaRQ(L+). 

This paper is organized as follows. In Section we highlight and solve numerical challenges that 
practitioners may face when implementing theoretical solutions for VaRa(T'’') in the homogeneous 
case Fi = ■ ■ ■ = Fd- Section presents the main concept underlying the Rearrangement Algorithm 
(RA) for computing Vail^(L+) and VaRQ(T“'“), levies criticism on its tuning parameters and 
investigates its empirical performance using various test cases. Section]^ then presents a conceptually 
and numerically improved version of the RA, which we call the Adaptive Rearrangement Algorithm 
(ARA), for calculating VaR^ .fL'*') and VaRa(L'''). Section [^concludes. Proofs of results stated in 
the main body are relegated to an appendix. 


2 Known optimal solutions in the homogeneous case and their 
tractability 

In order to assess the quality of general algorithms such as the RA, we need to know (at least some) 
optimal solutions with which we can compare such algorithms. Embrechts et al. [11[ Proposition 4] 
and Embrechts et al. Proposition 1] present formulas for obtaining VaRQ(T'’') in the homogeneous 
case. In this section, we address the corresponding numerical aspects and algorithmic improvements. 
We assume d > 3 throughout; for d = 2, Embrechts et al. EH Proposition 2] provide an explicit 
solution for computing VaRa(L'’') (under weak assumptions). 


2.1 Crude bounds for any VaRQ,(L+) 


The following lemma provides (crude) bounds for VaRQ(L+) which are useful for computing initial 
intervals (see Section 2.2) and conducting sanity checks. Note that we do not make any (moment or 
other) assumptions on the involved marginal loss distributions; in fact, they do not even have to be 
equal. Furthermore, the bounds do not depend on the underlying unknown copula. 


Lemma 2.1 (Crude bounds for VaRQ(T+)) 

Let Lj ^ Fj, j G {1, ... ,d}. For any a G (0,1), 


drain F- (a/d) < VaRa(T“'") < dmaxF- (l- - —V 

j J j \ d ■' 


( 1 ) 


where F- denotes the quantile function of Fj. 
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The bounds Q can be computed with the function crude_VaR_bounds () in the R package 
qrmtools. 

2.2 The dual bound approach for computing VaRQ(L+) 

This approach for computing VaRQ(T+) in the homogeneous case with margin(s) F is presented in 


Embrechts et al. 11, Proposition 4] and termed the dual bound approach in what follows; note that 
there is no corresponding algorithm for computing VaR ^,(L~^) with this approach. In the remaining 
part of this subsection we assume that T(0) = 0, F{x) < 1 for all x G [0, oo) and that F is absolutely 
continuous with ultimately decreasing density. Let 


D{s,t) = 


d 


s — dt. 


f'S—{d—l)t _ 

F{x)dx and 


Dis) = min D(s,t), 

t&[0,s/d] 


( 2 ) 


where F denotes the survival function of F, i.e., F{x) = 1 — F{x). In comparison to Embrechts et al. 

Proposition 4], the dual bound D here uses a compact interval for t (and thus min{-}) by our 
requirement E(0) = 0 and since lim^^^/^ H(s, t) = dF{s/d) by rHospital’s Rule. The procedure for 
computing VaRci(L''') according to Embrechts et al. 11 Proposition 4] can now be given as follows. 


Algorithm 2.2 (Computing VaRQ(L+) according to the dual bound approach) 

1) Specify initial intervals [s/,Su] and [tuty\. 

2) Inner root-finding in t: For each considered s G [s/, s„], compute D{s) by iterating over t G [ti,tu] 
until a t* is found for which h{s,t*) = 0, where 

h{s, t) := D{s, t) — {F{t) + {d — l)F{s — {d — l)f)). 

Then D{s) = F{t*) + {d - l)F{s - {d - 

3) Outer root-finding in s: Iterate Stepover s G [si, s„] until an s* is found for which D{s*) = 1—a. 
Then return s* = VaRQ(L'*“). 


Algorithm 2.2 is implemented in the function worst_VaR_hom(. . . , method="dual") in the R 
package qrmtools; the dual bound D is available via dual_bound() . It requires a one-dimensional 
numerical integration (unless F can be integrated explicitly) within two nested calls of a root-finding 
algorithm (uniroot () in R). Note that Algorithm |2.2| requires specification of the two initial intervals 
[si , Su] and [ti , and Embrechts et al. give no practical advice on how to choose them. 

First consider [ti,tu]. By our requirement E(0) = 0 one can choose ti = 0 (or the infimum of the 
support of F). For tu one would like to choose s/d; see the definition of D{s) in (|^. However, care 
has to be taken as h{s, s/d) = 0 for any s and thus the inner root-hnding procedure will directly 
stop when a root is found at tu = s/d. To take care of this, the inner root-hnding algorithm in 
worst_VaR_hom(. . . , method="dual") hxes f .upper, i.e., h{s,tu) for the considered s, to —/i(s,0) 
so that a root below tu = s/d can be detected; note that this is a (inelegant; for lack of a better 
method) adjustment in the function value, and not in the root-hnding interval [ti,tu]. 

Now consider [sj,Su], in particular, s;. According to Embrechts et al. 11 Proposition 4] si has to 
be chosen “sufficiently large”. If it is chosen too small, the inner root-hnding procedure in Step of 


Algorithm 2.2 will not be able to locate a root; see also the left-hand side of Figure There is 
currently no (good) solution known on how to automatically determine a sufficiently large si (note 
that worst_VaR_hom(. . . , method="dual") currently requires to be specihed by the user). 

Given si, one can then choose Su as the maximum of s; -|- 1 and the upper bound on VaRa(L''') as 
given in (Q, for example. 

On the theoretical side, the following proposition implies that if F is strictly convex, so is D{s, •) 
for fixed s. This shows the uniqueness of the minimum when computing D{s) as in (|^. A standard 
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result on convexity of objective value functions then implies that D{s) itself is also convex; see 
Rockafellar and Wets [21[ Proposition 2.22] and also the right-hand side of Figure]^ below. 

Proposition 2.3 (Properties of D{s,t) and D{s)) 

1) D{s) is decreasing. 

2) If F is convex, so is D(s,t). 


Example 2.4 (Auxiliary functions for the dual bound approach) 

As an example, consider d = 8 Par(0) risks with distribution function Fj{x) = 1 — (1 -f x)~^^ , a: > 0, 
6j >0. The left-hand side of Figure illustrates t i—)• h{s,t) for 6 = 2 and various s. Note that 
h{s, s/d) is indeed 0 and for (too) small s, h{s, t) does not have a root for t £ [0,s/d) as mentioned 
above. The right-hand side of Figure]^ shows the decreasing dual bound D{s) for various parameters 

6 . 




t 


S 


Figure 1 1 1— >■ h{s,t) for various s, d = 8 and F being Par(2) (left-hand side). The dual bound D{s) 
for d = 8 and F being Par(0) for various parameters 9 (right-hand side). 


2.3 Wang’s approach for computing VaRa(T+) 

Theory 

The approach mentioned in Embrechts et al. Proposition 1] is termed Wang’s approach here. 
It originates from Wang et al. |2^ Corollary 3.7] and, thus, strictly speaking, precedes the dual 
bound approach of Embrechts et al. ED Proposition 4]. It is conceptually simpler and numerically 
more stable than the dual bound approach, yet still not straightforward to apply. For notational 
simplicity, let us introduce 

ac = a + (d — l)c, 5c = 1 — c, (3) 

for c e [0, (1 — a)/d] (so that Oc G [a, 1 — (1 — ot)/d\ and 5c G [1 — (1 — a)/d, 1]) and 

I{c ):=-^— / F~{y)dy, c G (0, (1 - a)/d], 

with the assumption that F admits a density which is positive and decreasing on [f3, oo) for some 
P<F-{a). Then, for L~F, 


VaR„(L+) =dE[LlLG [F-(ac),E-(5c)]], a£[F{/3),l), 


( 4 ) 
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where c (typically depending on d, a) is the smallest number in (0, (1 — a)/d] such that 

I(c)>^F-(ae) + ^F-(6,). (5) 

In contrast to what is given in Embrechts et al. [^, note that (0, (1 — a)/d] has to exclude 0 since 
otherwise, for Par(0) margins with 0 G (0,1], c equals 0 and thus, erroneously, VaRQ,(L'’') = oo. If 
F is the distribution function of the Par(0) distribution, I is given by 


1(c) 


.6;^log(^) - 1, 


if 0 / 1, 
if 0 = 1. 


( 6 ) 


The conditional distribution function of L | L G [E (oc), F (be)] is given by 

^L|Le[F-(ac),F-(fec)](®) = ^ ^ [-^“(«c),E“(6c)]. Using this fact and a substitution, we 

obtain that, for a G [F{(3), 1), ^ becomes 


_ fF-ib.) 

VaR„(T+) = d / = d 

JF-(ac) 




be Ojc 


= dl{c). 


(7) 


Equation ([^ has the advantage of having the integration in 7(c) over a (at least theoretically) 
compact interval. Furthermore, finding the smallest c such that ^ holds also involves 7(c). A 
procedure thus only needs to know the quantile function F~ to compute VaRQ,(L'''). This leads to 
the following algorithm. 

Algorithm 2.5 (Computing VaRQ(L+) according to Wang’s approach) 

1) Specify an initial interval [ci,Cu] with 0 < c; < c„ < (1 — a)/d. 

2) Root-finding in c: Iterate over c G [c;,c^j] until a c* is found for which h{c*) = 0, where 

h{c) := 7(c) - (^E-(ae) + c G (0,(1 -«)/d]. (8) 


Then return {d — 1)E (oc*) -|- F {be*). 


This procedure is implemented in the function worst_VaR_hom(. . . , method="Wang") in the R 
package qrmtools with numerical integration via R’s integrate () for computing 7(c); the function 
worst_VaR_hom(. . . , method="Wang.Par") makes use of ([^. 

The following proposition shows that the root-finding problem in Step 2) of Algorithm 2.5 is 
well-defined in the case of Pareto margins for all 0 > 0 (including the infinite-mean case 
distributions under more restrictive assumptions, see bernardjiangwang2014 


for other 


Proposition 2.6 

Let F{x) = 1 — (1 -|- x)~^, 9 > 0, he the distribution function of the Par(0) distribution. Then h in 
Step 1^ of Algorithm |2.5| has a unique root on (0, (1 — a)/d), for all a G (0,1) and d > 2. 


Practice 

Let us now focus on the case of Par(0) margins (see worst_VaR_hom(. . . , method="Wang.Par")) 


in Algorithm 


2.5 


We first consider ci. I 
F at confidence 


and, in particular, how to choose the initial interval [c;,c 

satisfies 7(0) = ^~{y) dy = ESq,(L), i.e., 7(0) is the expected shortfall of L 

level a. If L has a finite first moment, then 1(0) is finite. Therefore, h{0) is finite (if E~(l) < oo) 
or —oo (if E~(l) = oo). Either way, one can take c; = 0. However, if L ~ E has an infinite first 
moment (see, e.g., Hofert and Wiithrich |17| or Chavez-Demoulin et al. for situations in which 
this can happen), then 7(0) = oo and E“(l) = oo, so h{0) is not defined; this happens, e.g., if E is 
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Par(0) with 0 e (0,1]. In such a case, we are forced to choose c/ e (0, (1 — a)/fi); see the following 
proposition for how this can be done theoretically. Concerning Cu, note that I’Hospital’s Rule implies 
that I{cu) = F~{1 — (1 — a)/d) and thus that /i((l — a)/d) = 0. We thus have a similar problem (a 
root at the upper endpoint of the initial interval) as for computing the dual bound. However, here 
we can construct a suitable < (1 — a)/d] see the following proposition. 


Proposition 2.7 (Computing ci,Cu for Par(6<) risks) 

Let F be the distribution function of a Par(0) distribution, 9 > 0. Then q and Cu in Algorithm 2.5 
can be chosen as 


Cl 


(l-e)(l-a) 

d ’ 

_1—a_ 

(d+i)^+d-F 

_1—a_ 

(d/(0-l)+l)e+d-l’ 


if 0G(O,1), 
if 0 = 1, 

if 0 G (1,oo). 


Cu — 


(l-a)(d-l+e) 
(d-l){2e+d) ’ 

1—a 

3(i/2-l’ 


if 0 / 1, 
if 0 = 1. 


In the following examples we briefly consider some experiments which have led to several numerical 
hurdles we had to overcome when implementing worst_VaR_hom(. . . , method="Wang.Par") ; see 
also the detailed vignette VaR_bounds for how to reproduce them (e.g.. Section 1.4). 

Example 2.8 {h, VaR ^,(L~*~) and VaRa(L+) for Wang’s approach with Par(0) risks) 

As an example, consider Par(0) risks and the confidence level a = 0.99. Figure illustrates the 
objective function h{c) as a function of c G (0, (1 — a)/d] (see ([^) for various 9 and d = 8 (left-hand 
side) and d = 100 (right-hand side). Non-positive values h{c) have been omitted so that the y-axis 
could be given in log-scale; note how steep h can be (we have evaluated h at 2^^ -|- 1 points equally 
spaced between 0 and (1 — a)/d), especially for small 9 (and large d). For an even larger number of 
evaluation points c, one sees that the first and last positive values of h indeed reach down towards 
0. Overall, the objective function can be evaluated without numerical problems on (0, (1 — aj/d] for 
our chosen 9 and d. 

Figure [^displays VaR„ (L~*~) and VaRQ,(L+) as a function in 1 — a for various 9 and d = 8 
(left-hand side) and d = 100 (right-hand side). 




C C 

Figure 2 Objective function h{c) for a = 0.99, F being Par(0), d = 8 (left-hand side) and d = 100 
(right-hand side). 


For obtaining numerically reliable results (over these wide ranges of parameters; indeed we 
tested much higher dimensions as well), one has to be careful when computing the root of h for 
c G (0, (1 — a)/d). First, choosing a smaller root-finding tolerance is crucial. Figurebelow (see also 
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Figure 3 VaR^ /L"*") and VaRQ(L''') as functions of 1 — a for F being Par(0), d = 8 (left-hand side) 
and d = 100 (right-hand side). 


Example |2.9[ ) shows what silently happens if this is not considered (our procedure chooses MATLAB’s 
default 2.2204 • 10“^® instead of the much larger unirootO default .Machine$double.eps''0.25). 
Second, it turned out to be required to adjust the theoretically valid initial interval described in 
Proposition |2 . 7| further in order to guarantee that h is numerically of opposite sign at the interval 
end points. In particular, worst_VaR_hom(. . ., method="Wang.Par") chooses c;/2 as lower end 
point (with c; as in Proposition 2.7) in case 6^1. 

These problems are described in detail in Section 1.4 of the vignette VaR_bounds, where we also 
show that transforming the auxiliary function h to a root-finding problem on (l,oo) as described in 
the proof of Proposition |2.7| does not require a smaller root-finding tolerance but also an extended 
initial interval and, furthermore, faces a cancellation problem (which can be solved, though); 
see also the left-hand side of Figure where we compare this approach to worst_VaR_hom(. . . , 
method="Wang.Par") after fixing these numerical issues. 

In short, one should be very careful when implementing supposedly “explicit solutions” for 
computing VaR^(L+) or VaRa(L^) in the homogeneous case with Par(0) (and most likely also 
other) margins. 


Example 2.9 (Comparison of the approaches for Par(0) risks) 

Again let us consider Par(0) risks and the confidence level a = 0.99. Figure compares Wang’s 
approach (using numerical integration; see worst_VaR_hom(. . . , method="Wang")), Wang’s ap¬ 
proach (with an analytical formula for the integral J(c) but uniroot ()’s default tolerance; see 
the vignette VaR_bounds), Wang’s approach (with an analytical formula for the integral /(c) and 
auxiliary function h transformed to (l,oo); see the vignette VaR_bounds), Wang’s approach (with 
analytical formula for the integral /(c), smaller unirootO tolerance and adjusted initial interval; 
see worst_VaR_hom(. . . , method= "Wang.Par")), and the lower and upper bounds obtained from 
the RA (with absolute tolerance 0); see Section]^ and RA(). All of the results are divided by the 
values obtained from the dual bound approach to facilitate comparison. The two plots (for d = 8 
and d = 100, respectively) show that comparable results are obtained by the different approaches 
and why it is important to use a smaller tolerance for unirootO in Wang’s approach. 


Let us again stress how important the initial interval [c;, Cu] is. One could be tempted to simply 
choose c„ = (1 — a)// and force the auxiliary function h to be of opposite sign at c^, e.g., by setting 






























Standardized (by dual method) VaRo.gg for d = 8 Par(9) margins 
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Figure 4 Comparisons of Wang’s approach (using numerical integration; see worst_VaR_hom(. . . , 
method="Waiig") ), Wang’s approach (with an analytical formula for the integral /(c) but 
unirootO’s default tolerance; see the vignette VaR_bouiids), Wang’s approach (with an 
analytical formula for the integral /(c) and auxiliary function h transformed to (l,oo); 
see the vignette VaR_bounds), Wang’s approach (with analytical formula for the integral 
/(c), smaller unirootO tolerance and adjusted initial interval; see worst_VaR_hom(. . . , 
method="Wang.Par")), and the lower and upper bounds obtained from the RA; all of the 
results are divided by the values obtained from the dual bound approach to facilitate 
comparison. The left-hand side shows the case d = 8, the right-hand side d = 100. 
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h{cu) to .Machine$double .xmin, a positive but small number. Figure shows graphs similar to 
the left-hand sides of Figures (but VaRQ,(L+) only) and (standardized with respect to the upper 
bound obtained from the RA). In particular, VaRQ.(L'’') is not monotone in a anymore (see the 
left-hand side of Figure |5[) and the computed VaRQ(L''') values are not correct anymore (see the 
right-hand side of FigurelM. 
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Figure 5 Figures corresponding to the left-hand side of Figures (VaRQ(F“'“) only) and [^(stan¬ 
dardized with respect to the upper bound obtained from the RA) but for h{{l — a)/d) 
adjusted to .Machine$double.xmin. 


After carefully considering all the numerical issues, we can now look at VaR^ ,(L^l and VaRQ,(L''') 
from a different perspective. The right-hand side of Figure]^ shows VaR ^(L^) and VaRQ(L'*“) as 
functions in the dimension d. The linearity of VaRQ(L'’') in the log-log scale suggests that VaRQ(T'’') 
is actually a power function in d. To the best of our knowledge, this is not known (nor theoretically 
justified) yet. 


3 The Rearrangement Algorithm 

We now consider the RA for computing VaR q(L'*~) and VaRa(L''') in the inhomogeneous case; as 
before, we mainly focus on VaRQ(L+) here. 

3.1 About the algorithm 

One of the early works on finding bounds for VaRci(L'*“) (including a proof of their sharpness) can 
be found in Makarov |18] , who provides bounds on the distribution function of the sum of two 
random variables and bounds on VaRQ,(L'*') thus follow. Later on, Firpo and Ridder |14] prove 
these results using copula theory, introducing dependence structures into the above framework 
and extend Makarov’s results to include an arbitrary increasing continuous aggregation function 
(they do not prove the sharpness of the bounds, though). Williamson and Downs develop new 
methods for calculating convolutions and dependency bounds for the distributions of functions of 
random variables; they use lower and upper approximations to the desired distribution, containing 
the representation error, and provide bounds on the errors. 
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Figure 6 A comparison of VaRQ(A'*') computed with worst_VaR_hom(. . . , method="Wang.Par") 
(solid line) and with the approach based on transforming the auxiliary function h to a, 


root-finding problem on (1, oo) (dashed line) as described in the proof of Proposition 2.7 


(left-hand side). VaR ,^,(L'*~l (dashed line) and VaRQ(L'’“) (solid line) as functions in d on 
log-log scale (right-hand side). 


Denuit et al. extend the above two-dimensional frameworks and show how to compute bounds 
on the distribution function of L"*" = Li for the d >3 case. In a similar work, Cossette et al. 

further develop results about {L\,L 2 ) by assuming additional information about the correlation 
structure of {Li,L 2 ). They further extend their results to the general multivariate case and propose 
bounds for continuous and componentwise monotone functions in Li,..., assuming that the 
only available information on Lj is its distribution function Fj, j € {1,... ,d}. By relaxing some of 
the continuity assumptions with respect to the aggregation function, Embrechts et al. provide a 
generalization of these results using copula theory. In the latter paper, it is shown that without any 
prior information on the dependence structure, only bounds for the distribution function of the sum 
of the risks can be found and the problem of the sharpness of these bounds remained open when 
d>2. 


Embrechts and Puccetti provide better bounds on in the homogeneous case Fi = ■ ■ ■ = 
Fd = F, see Section 2.2 based on the duality result of Riischendorf 23 for a continuous distribution 
function F. The assumption of equal margins is rather restrictive, especially for large d. To address 
this problem, Embrechts and Puccetti extend the dual bound approach to general portfolios and 
describe a numerical procedure to compute a lower bound for Fi^+ for an inhomogeneous portfolio 
of risks. The shortcoming of this method is that it requires the application of a global optimization 
algorithm for which there is no guarantee of convergence to a global optimum in a reasonable 
and predictable amount of time; more importantly the quality of performance of many of these 
optimization procedures is not well understood and in general the performance of such algorithms 
deteriorates as d increases. For these reasons the application of this method for d > 50 becomes 
intractable in some cases. 


The above problems can be studied in the context of completely mixable matrices, i.e., matrices 
such that there exists a collection of permutations acting on the columns which result in a constant 
row sum. This concept of complete mixability is introduced and discussed in Wang and Wang 
(^. If a matrix is not completely mixable, determining the smallest maximal and largest minimal 
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row sums is of interest to find VaR^ ,(L^) and VaRQ(L+), respectively (or at least approximations 
of such). These minimax and maximin problems, respectively, are in turn connected to discrete 
approximations of the marginal quantile functions. Haus |15| points out that such approaches for 


computing VaR^(L+) and VaRQ,(L+) are related to the multidimensional bottleneck assignment 
problem and complete mixability is in general AfR-complete. As a result, convergence to an optimal 
solution in polynomial time is not guaranteed. 

If we discretize the tail of each of the d risk factors using N points, the underlying space over 
which the above problems are analyzed becomes an N x d matrix and enumeration of the total 
number of possible matrices arising from permuting all but one column becomes intractable in 
applications as there are possible matrices. One can easily observe that for a portfolio of 

only 10 positions and N = 20, the total number of such matrices is (20!)® £ 0(10^®^). Note that 
the choices of both = 20 and d = 10 are extremely conservative and for illustrative purposes only; 
in practice N can easily be as large as 1000 to 100,000 and many portfolios consist of at least 20 to 
40 instruments; larger financial institutions can have portfolios with several thousand instruments. 

As noted earlier, the complexity of the optimization procedures required for finding dual bounds, 
given arbitrary marginal distributions along with a high run time were the main drawbacks of 
using the dual bounds approach for obtaining a lower and upper bound for VaRo(L''“). As a 
result, Puccetti and Riischendorf propose the Rearrangement Algorithm (RA) to tackle these 
issues. The initial idea underlying the RA and the numerical approximation introduced in Puccetti 
and Riischendorf 
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for calculating VaR„ y(L~*~) and VaRa(T“'") is due to Riischendorf 
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and 

Riischendorf j^, respectively. The higher accuracy of the RA (theoretically still an open question) 
along with its simple implementation compared to the previous methods makes the RA an attractive 
alternative for obtaining VaR .,^,(L'*~) and VaRQ(T'’') when (only) the marginal loss distributions 
Fi,... ,Fci are known. In the following subsections we look at how the RA works and analyze its 
performance using various test cases. 


3.2 How the Rearrangement Algorithm works 

The RA can be applied to approximate the best Value-at-Risk VaRg lL'*') or the worst Value-at-Risk 
VaRQ(L'’') for any set of marginals Fj, j € {1, ... ,d}. In what follows our focus is mainly on 
VaRQ(L'''); our implementation RA() in the R package qrmtools also addresses VaR„(L+). To 
understand the algorithm, note that two columns a, 6 £ IR'^ are called oppositely ordered if for 
all i,j £ {1,..., N} we have (oj — aj){bi — bj) < 0. Given a number N of discretization points 


below), the RA constructs two [N, (i)-matrices, denoted by A“ and X ; the hrst matrix aims at 
constructing an approximation of VaRo(L'’“) from below, the second matrix is used to construct 
an approximation of VaRa(T^) from above. Separately for each of these matrices, the RA iterates 
over its columns and permutes each of them so that it is oppositely ordered to the sum of all other 
columns. This is repeated until the minimal row sum 


of the marginal quantile functions F^ ,... ,F^ above a (see Steps 2.1) and 3.1) of Algorithm 


3.1 


s(A) = min > : 

l<t<N 


(for X = {xij) being one of the said (N, d)-matrices) changes by less than a given (convergence) 
tolerance e > 0. The RA for VaRQ,(L^) thus aims at solving the maximin problem. The intuition 
behind this is to minimize the variance of the conditional distribution of > F£^{a) to 

concentrate more of the 1 — a mass of Fi+ in its tail. This pushes VaRa(L+) further up. As 
Embrechts et al. state, one then typically ends up with two matrices whose minimal row sums 
are close to each other and roughly equal to VaRQ,(L'’'). Note that if one iteration over all columns 
of one of the matrices does not lead to any change in that matrix, then each column of the matrix 
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is oppositely ordered to the sum of all others and thus there is also no change in the minimal row 
sum (but the converse is not necessarily true, see below). 

The version of the RA given below contains slightly more information than in Embrechts et ah 
111]; e.g., how inhnite quantiles are dealt with. For more features of the actual implementation 


which is an improved version of the one given below (see Section 3.3), see RA() and the underlying 
workhorse rearrange (). This includes, e.g., a parameter max.ra determining the maximal number 
of columns to rearrange or the choice e = (abstol=)NULL to iterate until each column is oppositely 
ordered to the sum of all others. The latter is typically (by far) not implied by e = 0, but does 
not bring better accuracy (see, e.g., the application discussed in the vignette VaR_bounds) and is 
typically very time-consuming (hence the introduction of max. ra). 


Algorithm 3.1 (RA for computing VaRQ(T^)) 

1) Fix a confidence level a G (0,1), marginal quantile functions Ff 
the desired (absolute) convergence tolerance e > 0. 

2) Compute the lower bound: 


F, , an integer N € N and 


2.1) Define the matrix = (xfj) for xfj = F~ {a + i G {1,..., Ai}, j G {1,... ,d}. 

2.2) Permute randomly the elements in each column of 

2.3) For 1 < j < d, permute the j-th column of the matrix A" so that it becomes oppositely 
ordered to the sum of all other columns. Call the resulting matrix Y°‘. 


2.4) Repeat Step 2.3) until s(Y°^) — s(A“) < s, then set = s(y“). 

3) Compute the upper bound: 

3.1) Dehne the matrix = (xfj) for = F~ (a -|- ^^^^), i G {1,..., N}, j G {1,..., d}. If 

(for i = N and) for any j G {1,... ,d}, F~{1) = oo, adjust it to F~ {a -\- _ 

3.2) Permute randomly the elements in each column of . 

3.3) For 1 < J < d, iteratively rearrange the j-th column of the matrix X"^ so that it becomes 
oppositely ordered to the sum of all other columns. Call the resulting matrix Y . 


3.4) Repeat Step 3.3) until s{Y ) — s{X ) < e, then set sat = s(y ). 

4) Return (s^r, sn)- 

As mentioned before, the main feature of the RA is to iterate over all columns and oppositely 
order each of them with respect to the sum of all others (see Steps 2.3)| and 3.3) |. This procedure 
aims at reducing the variance of the row sums with each rearrangement. Note that it does not 
necessarily reach an optimal solution of the maximin problem (see, e.g., Haus |15[ Lemma 6] for a 
counter-example) and thus the convergence \sn — —>■ 0 is not guaranteed. In order to reduce the 

possibility of this happening in practice, the randomization of the initial input in Steps 2.2)|and 3.2) 


has been put in place; however, see our study in Section 2.3 of the vignette VaR_bounds concerning 
the possible rearranged output matrices and the influence of the underlying sorting algorithm on 
the outcome. 


3.3 Conceptual and numerical improvements 

Some words of warning are in order here. Besides the confidence level a and the marginal quantile 
functions F ^,..., F^, RA relies on two sources of input, namely A G N and e > 0, for which 
Embrechts et al. jllj do not provide guidance on reasonable defaults. Concerning N, it obviously 
needs to be “sufficiently large”, but a practitioner is left alone with such a choice. Another issue is 
the use of the absolute tolerance e in the algorithm. There are two problems. The first problem is 
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that it is more natural to use a relative instead of an absolute tolerance in this context. Without 
(roughly) knowing the minimal row sum in Steps 2.4) and 3.4) a pre-specihed absolute tolerance 
does not guarantee that the change in the minimal row sum from to is of the right order 
(and such order depends at least on d and the chosen quantile functions). If e is chosen too large, 
the computed bounds and sn would carry too much uncertainty, whereas if it is too small, an 
unnecessarily long run time results; the latter seems to be the case for Embrechts et al. [IT| Table 3], 
where the chosen e = 0.1 is roughly 0.000004% of the computed VaRo, 99 (T+). 

The second problem is that the absolute tolerance e is only used for checking individual “conver¬ 
gence” of Sjv and of sat. It does not guarantee that and sn are sufficiently close to obtain a 
reasonable approximation to VaRQ,(L'’'). We are aware of the theoretical hurdles underlying the 
algorithm which are still open questions at this point (e.g., the probability of convergence of 
and sn to VaRQ(T"'') or that VaRo(L+) <sn for sufficiently large N), but from a computational 
point of view one should still check that are close to each other. Also, the algorithm 

should return convergence and other useful information, e.g., the relative rearrangement range 
|(sAr — Sj^)/sn\, the actual individual absolute tolerances reached when computing and sn, 
the number of column rearrangements used, logical variables indicating whether the individual 
absolute tolerances have been reached, the number of column rearrangements for and sat, the 
row sums computed after each column rearrangement, the constructed input matrices A“, X and 
the corresponding rearranged, final matrices Y^, Y ; see RA() in the R package qrmtools for such 
information. 

Another suboptimal design of the RA is to iterate over all d columns before checking the termination 
conditions; see Steps 2.3) and 3.3) of Algorithm 3.1 Our underlying workhorse rearrange () keeps 
track of the column rearrangements of the last d considered columns and can thus terminate after 
rearranging any column (not only the last one); see also Algorithm 4.1 below. This saves run time 
(despite the “tracking” overhead). We advise the interested reader to have a look at the source code 
of rearrange () for further numerical and run-time improvements (fast accessing of columns via 
lists; avoiding having to compute the row sums over all but the current column; an extended tracing 
feature), some of which are mentioned in the vignette VaR_bounds. 


3.4 Empirical performance under various scenarios 

In order to empirically investigate the performance of the RA, we consider two studies, each of 

which addresses four cases; we thus consider eight scenarios. As studies, we consider the following: 

Study 1: N £ {2\ 2^ ..., 2^'^} and d = 20; 

Study 2: Ai = 2® = 256 and d E {2^, 2®,..., 2^°}. 

These choices allow us to investigate the impact of the upper tail discretization parameter N (in 

As cases, we consider the following different marginal tail behaviors based on the Pareto distribution 

function Fj{x) = 1 — {1 + 

Case HH : 9i,... ,9d form an equidistant sequence from 0.6 to 0.4; this case represents a portfolio 
with all marginal loss distributions being heavy-tailed (slightly increasing in heavy- 
tailedness). 

Case LH : 0i,..., form an equidistant sequence from 1.5 to 0.5; this case represents a portfolio 
with different marginal tail behaviors ranging from comparably light-tailed to very 
heavy-tailed distributions. 

Case LL : 9i,... ,9d form an equidistant sequence from 1.6 to 1.4; this case represents a portfolio 
with all marginal loss distributions being comparably light-tailed. 


Study 1) and the impact of the number of risk factors d (in Study 2) on the performance of the RA. 
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Case LHi : 6i,..., 6d-i are chosen as in Case LL and 0d = 0.5; this case represents a portfolio all 
marginal loss distributions being light-tailed except the last. 

To keep the studies tractable, we focus on the confidence level a = 0.99 and the absolute convergence 
tolerance e = 0 in all scenarios. Furthermore, we consider B = 100 replicated simulation runs in 
order to provide empirical 95% confidence intervals for the estimated quantities; note that some of 
them are so tight that they are barely visible in the figures presented below. The B replications only 
differ due to different permutations of the columns in Steps [2^^ and [33)] of Algorithm |3.1[ everything 
else is deterministic; this allows us to study the effect of these (initial) randomization steps on the 
(convergence) results of the RA. Concerning the hardware used, all results were produced on an 
AMD 3.2 GHz Phenom II X4 955 processor with 8 GB RAM. 

Results of Study 1 {N running, d fixed) 

The simulation results for [Study l| can be summarized as follows: 

■ As can be seen in Figure the means over all B computed Sjy sn converge as N increases. 

■ Figure [^indicates that as N increases, so does the mean elapsed time (as to be expected). Overall, 
run time does not drastically depend on the case for our choices of Pareto margins, which is a 
good feature. 

■ Figure 1^ shows that the maximum number of column rearrangements rarely exceeds lOd as N 
increases; this will be used as a default maximum number of column rearrangements required in 
the ARA. 

■ Finally, Figure [^indicates that the rate at which the number of oppositely ordered columns 
(based on the hnal rearranged T" and Y°^) decreases depends on the characteristics of the 
marginal distributions involved, i.e., the input matrix X. The number of oppositely ordered 
columns seems particularly small (for large N) in Case LL, where essentially only the last column 
is oppositely ordered to the sum of all others. 


- Mean(upper bound) 

— Bootstrapped 95% CIs 
- Mean (lower bound) 

s ■ 

- - Bootslrap^d 95% CIs 
- Mean (lower bound) 

1200 

— Bootstrapped 95% CIs 
- Mean (lower bound) 

15800 



180000 18 


1160 1180 


1 - 

S 

'r 

_ 


r 


L 



_ r 

f Y 


f 5- 


f 1- 





I - 


5000 



X o - 

X - - 


I o 

_J O _ 


_l 

-1 S 

- Mean (upper bound) 

- - Bootstrapped 95% CIs 


o 


o 


o 

- - Bootstrapped 95% CIs 


Figure 7 Study I: VaRo .99 bounds Spf and sn for the Cases HH, LH, LL and LHi (from left to 
right). 


Results of Study 2 {N fixed, d running) 


Figures 11-14 show the performance of the RA in Study 2[ here we are interested in analyzing the 
impact of the number of risk factors d on portfolios which exhibit different marginal tail behaviors. 
The simulation results can be summarized as follows: 


■ Figureshows that the means over all computed Sjy and sjv diverge from one another; especially 
when more marginal distributions are heavy-tailed. This is due to the fact that we have kept N 
the same for all cases in|Stud^(2 
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Similar to what we have seen in [Study l] Figure 12 indicates that the mean run time of the 
RA increases as the number of risk factors increases; Case LL in Study 2 has the least run time 
on average as 9i,... ,9d form an equidistant sequence from 1.6 to 1.4 which results in smaller 
elements of the input matrix X with a smaller range of entries compared to the other cases. 


■ As in |StudyT!| the mean number of rearranged columns is typically below lOd; see Figure [T^ 

■ The number of oppositely ordered columns (of y“ and y°‘), see Figure [I^ increases as d increases; 
note that this finding does not contradict what we have seen in Study 1 as we have kept N the 
same here. 






Figure 11 Study 2: VaRo .99 bounds and sn for the Cases HH, LH, LL and LHi (from left to 
right). 




Figure 12 Study 2: Run times (in s) for the Cases HH, LH, LL and LHi (from left to right). 




Figure 13 Study 2: Number of rearranged columns for the Cases HH, LH, LL and LHi (from left 
to right). 
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Figure 14 Study 2: Number of oppositely ordered columns (of Yf‘ and Y°‘) for the Cases HH, LH, 
LL and LHi (from left to right). 


4 The Adaptive Rearrangement Algorithm 

4.1 How the Adaptive Rearrangement Algorithm works 

In this section we present an adaptive version of the RA, termed Adaptive Rearrangement Algorithm 
(ARA). This algorithm for computing the bounds sjy for VaR^ dL^l or VaRQ,(L''') (as 

before, we focus on the latter) provides an algorithmically improved version of the RA, has more 
meaningful tuning parameters and adaptively chooses the number of discretization points N. The 
ARA is implemented in the R package qrmtools, see the function ARA(). Similar to our RA() 
implementation, ARA() also relies on the workhorse rearrange () and retnrns mnch more information 
(and can also compute VaR „(L~*~ll, but the essential part of the algorithm is given as follows. 

Algorithm 4.1 (ARA for computing VaRQ(L+)) 

1) Fix a confidence level a G (0,1), marginal quantile functions an integer vector 

K G N^, I G N, (containing the numbers of discretization points which are adaptively used), a 
bivariate vector of relative convergence tolerances e = (ei,e 2 ) (containing the individnal relative 
tolerance ei and the joint relative tolerance 62 ; see below) and the maximal nnmber of iterations 
used for each k & K. 

2) For IV = 2^, A: G AT, do: 

2 . 1 ) Compute the lower bound: 

2.1.1) Define the matrix = (xg) for xfj = F~ (a + i G {1,... ,N}, j G 

{l,...,d}. 

2.1.2) Permute randomly the elements in each column of A". 

2.1.3) Iteratively rearrange the jth column (j G {1, 2,..., d, 1,2,..., d,... }) of X°‘ op¬ 
positely to the sum of all other columns until the maximal number of column 
rearrangements has been reached or until 


s(i:“) - s(A“) 


< ei, 


(9) 


where denotes the matrix of quantiles after the jth colnmn has been rearranged 
and A“ denotes the same matrix d steps earlier, i.e., after the jth column has been 
rearranged the last time. 

2.1.4) Set sj^ = s(i:“). 

2.2) Compute the upper bound: 
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2.2.1) Define the matrix = (xfj) for xfj = F- {a + i G {1,...,A^}, j G 

{l,...,h}. If (for i = N and) for any j G {l,...,(i}, .^^”(1) = oo, adjust it to 

2.2.2) Permute randomly the elements in each column of . 

2.2.3) Iteratively rearrange the jth column (j G {1, 2,..., d, 1,2,..., d,... }) of X°‘ op¬ 
positely to the sum of all other columns until the maximal number of column 
rearrangements has been reached or until 


s(y") - s(x“) 
s(X“) 


< ei, 


( 10 ) 


where denotes the matrix of quantiles after the jth column has been rearranged 
and X denotes the same matrix d steps earlier, i.e., after the jth column has been 
rearranged the last time. 

2.2.4) Set SN = s(F“). 

2.3) Determine convergence based on both the individual and the joint relative convergence 
tolerances: 


If Q and (10) hold, and if 


SN — Sj^ 


SN 


< £ 2 , break. 


3) Return {sj^, sn)- 


Concerning the choices of tuning parameters in Algorithm 4.1 note that if = (A: = log 2 N), 


so if we have a single number of discretization points, then the ARA reduces to the RA but uses 
the more meaningful relative instead of absolute tolerances and not only checks what we termed 
individual (convergence) tolerance, i.e., the tolerance ei for checking “convergence” of and of sat 
individually, but also the jomt (convergence tolerance), i.e., the relative tolerance 62 between and 
sjv; furthermore, termination is checked after each column rearrangement (which is also done by our 
implementation RA() but not the version given in Embrechts et al. (^). As our simulation studies 


in Section 3.4 suggest, useful (conservative) defaults for K and the maximal number of column 
rearrangements are K = {8,9,..., 19) and lOd, respectively. Given the high model uncertainty and 
the (often) rather large values of YaRa{L~^) (especially in heavy-tailed test cases), a conservative 
choice for the relative tolerance e may be e = ( 0 , 0 . 01 ); obviously, all these values can be freely 
chosen in the actual implementation of ARA (). 


4.2 Empirical performance under various scenarios 

In this section, we empirically investigate the performance of the ARA. To this end, we consider 


d G {20,100} risk factors, paired with the Cases HH, LH, LL, LHi as in Section 3.4 The considered 


relative joint tolerances are 0.5%, 1% and 2% and we investigate the results for the individual 
relative tolerances 0% and 0.1%. Therefore the performance of ARA is investigated in 48 different 
scenarios. As before, the results shown are based on B = 100 simulations and we investigate the 
VaRo. 99 (T'’') bounds and sn, the N used on the final column rearrangement of ARA (i.e., the 
N = 2^, k £ K for which the algorithm terminates), the run time (in s), the number of column 
rearrangements (measured for the N used on termination of the algorithm) and the number of 
oppositely ordered columns after termination of the algorithm. 
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Results 


We first consider the results for the individual relative tolerance chosen as ei = 0. 


Our findings (see Figures 15 18 and computed numbers; the latter are omitted) indicate that: 
Although for both d = 20 and d = 100 the length of the confidence intervals for VaRo, 99 (-h'’') can 
be checked to be increasing as the joint relative tolerance 82 gets larger, the mean and lower 
and upper confidence bounds remain fairly close to each other. More importantly for a fixed 
individual relative tolerance, as 82 increases, we do not observe a drastic shift in both lower and 
upper bounds for the mean across different examples. 

An important observation about the N in the ARA is that in virtually all examples, the 95% 
confidence interval remains the same; this fact can be leveraged in practice for portfolios which 
exhibit the same marginal tail behavior to reduce the run time of the ARA. 

Across all of the 24 scenarios, doubling the joint relative tolerance reduces the run time (measured 
in s) more than 50%. 

The number of column rearrangements for the N used remains below lOd. 


■ Finally, as Figures 15 -18 reveal, randomizing the input matrix X has a minimal impact on 
various outputs of the ARA. However, this randomization has an interesting effect on the run 
time in that it seems to avoid the worst case in which sorting a lot of numbers is required when 
oppositely ordering the columns causing the algorithm to take quite a bit longer. This behavior 
(and dependence of run time on the order of the columns as well; typically, convergence was 
faster with heavier-tailed columns last) was clearly visible during our testing phase and leads to 
the following open research question: How can the rows and columns of the initial matrices be 
sorted such that the run time of the ARA is minimal? 

The effect of a different choice of 81 can be seen from Figures |20f[^ Overall, they are very similar, 
note however Figure 25 (based on the 20 constituents of the SMI from 2011-09-12 to 2012-03-28) 
concerning a too large choice of 81 . As we can see, both bounds can be fairly close (according to 
82 ) but the individual “convergence” did not take place yet; hence our default ei = 0 of RA() and 
ARA(). 
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Figure 15 Boxplots of the lower and upper VaRo. 99 (T+) bounds Sjv (left-hand side) and sjy (right- 
hand side) computed with the ARA for itol=0 based on B = 100 replications. 


5 Conclusion 

This paper presents two major contributions to the computation of the worst Value-at-Risk for a 
sum of losses with given marginals in the context of Quantitative Risk Management. 
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HHLHLLLH, HHLHLLLH, HHLHLLLH, HHLHLLLH, HHLHLLLH, HHLHLLLH, HHLHLLLH, HHLHLLLH, HHLHLLLH, HHLHLLLH, HHLHLLLH, HHLHLLLH, 

E2 = 0.005 E2 = 0.01 82 = 0.02 £2 = 0-005 82=0.01 E2 = 0.02 £2 = 0.005 82 = 0.01 82 = 0.02 82 = 0.005 £2 = 0.01 £2 = 0.02 

d = 20 d = 20 c/ = 20 d = 100 d=100 d = 100 d = 20 d = 20 d = 20 d=100 d = 100 d = 100 


Figure 16 Boxplots of the actual N = 2^ used for computing the lower and upper VaRo. 99 (-h'^) 
bounds Sjv (left-hand side) and sn (right-hand side) with the ARA for itol=0 based on 
B = 100 replications. 
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Figure 17 Boxplots of the run time (in s) for computing the lower and upper VaRo. 99 (-h^) bounds 
(left-hand side) and sw (right-hand side) with the ARA for itol=0 based on R = 100 
replications. 






Improved Algorithms for Computing Worst VaR 
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Figure 18 Boxplots of the number of column rearrangements (measured for the N used) for com¬ 
puting the lower and upper VaRo. 99 (i''') bounds sjv (left-hand side) and sn (right-hand 
side) with the ARA for itol=0 based on B = 100 replications. 
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Figure 19 Boxplots of the number of oppositely ordered columns for computing the lower and upper 
VaRo. 99 (A+) bounds (left-hand side) and sn (right-hand side) with the ARA for 
itol=0 based on B = 100 replications. 
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Figure 20 Boxplots of the lower and upper VaRo. 99 (A+) bounds Sjy (left-hand side) and sn (right- 
hand side) computed with the ARA for itol=0.001 based on B = 100 replications. 
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Figure 21 Boxplots of the actual N = 2^ used for computing the lower and upper VaRo. 99 (-h+) 
bounds Sjv (left-hand side) and sn (right-hand side) with the ARA for ito 1=0.001 based 
on B = 100 replications. 
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Figure 22 Boxplots of the run time (in s) for computing the lower and upper VaRo. 99 (-h^) bounds 
S]\f (left-hand side) and sn (right-hand side) with the ARA for itol=0.001 based on 
B = 100 replications. 
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Figure 23 Boxplots of the number of column rearrangements (measured for the N used) for com¬ 
puting the lower and upper VaRo. 99 (i''') bounds (left-hand side) and sw (right-hand 
side) with the ARA for itol=0.001 based on B = 100 replications. 
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Figure 24 Boxplots of the number of oppositely ordered columns for computing the lower and upper 
VaRo. 99 (A+) bounds (left-hand side) and sm (right-hand side) with the ARA for 
itol=0.001 based on B = 100 replications. 
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Relative tolerance tol of rearrange() 

Figure 25 Computed VaRo, 99 (-^^’'") bounds 

and sn (with rearrange()) depend¬ 
ing on the chosen relative tolerance 
£i (the cross “-I-” indicating the val¬ 
ues for El = NULL) for an application 
to SMI constituents data from 2011- 
09-12 to 2012-03-28; see the vignette 
VaR bounds for more details. 



Figure 26 Relative speed-up (in %) of the ac¬ 
tual implementation of ARA () over a 
basic version as given in the vignette 
VaR bounds. 


First, we considered the homogeneous case (i.e., all margins being equal) and addressed the 
dual bound approach based on Embrechts et al. (m Proposition 4] and Wang’s approach based 
on Embrechts et al. Proposition 1] for computing worst Value-at-Risk. Although both of these 
approaches are mathematically “explicit”, care has to be exercised when computing worst Value-at- 
Risk bounds with these algorithms in practice. We identified and overcame several numerical and 
computational hurdles in their implementation and addressed them using the R package qrmtools 
including the vignette VaR_bounds. We covered several numerical steps such as how to compute 
initial intervals for the root-finding procedure involved; a particular example which highlights the 
numerical challenges when computing worst Value-at-Risk in general is the case of equal Pareto 
margins (for which we also showed uniqueness of the root even in the infinite-mean case). 

Overall, the reader should keep in mind that there is a substantial difference between implementing 
a specific model (say, with Par(2) margins) where initial intervals can be guessed or chosen “suffi¬ 
ciently large” and the proper implementation of a result such as Embrechts et al. [IT| Proposition 4] 
in the form of a black-box algorithm; see the source code of qrmtools for more details and the work 
required to go in this direction. 

Second, we considered the general, i.e., inhomogeneous case. We first investigated the Rear¬ 
rangement Algorithm presented by Embrechts et al. [11] . Due to its simplicity, this algorithm 
by now has been widely adopted by the industry (see also https://sites.google.com/site/ 
rearrangementalgorithm/). Nevertheless, the original algorithm leaves questions unanswered 
concerning the concrete choice of two of its tuning parameters. These parameters were shown to 
have a substantial impact on the algorithm’s performance and thus need to be chosen with care. We 
therefore presented an improved version of the Rearrangement Algorithm termed Adaptive Rear¬ 
rangement Algorithm. The latter improves the former in that it addresses the aforementioned two 
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tuning parameters and improves on the underlying algorithmic design. The number of discretization 
points is chosen automatically in an adaptive way (hence the name of the algorithm). The absolute 
convergence tolerance is replaced by two relative convergence tolerances. Since they are relative 
tolerances, their choice is much more intuitive. Furthermore, the first relative tolerance is used to 
determine the individual convergence of each of the lower bound Sj^ and the upper bound sjsf for 
worst Value-at-Risk and the second relative tolerance is used to control how far apart and sn 
are; the original version of the Rearrangement Algorithm does not allow for such a control. The 
Adaptive Rearrangement algorithm has been implemented in the R package qrmtools, together 
with conservative defaults. The implementation contains several other improvements as well (e.g., 
fast accessing of columns via lists; avoiding having to compute the row sums over all but the current 
column; an extended tracing feature). 

There are still several interesting questions left to be investigated. First of all, as for the 
Rearrangement Algorithm, the theoretical convergence properties of the Adaptive Rearrangement 
Algorithm remain an open problem. Also, as mentioned above, it remains unclear how the rows and 
columns of the input matrices for the RA or ARA can be set up in an initialization such that the 
run time is minimal. Another interesting question is by how much we can reduce run time when 
using the rearranged matrix from case N = to construct the matrix for the case N = 2^ (which 
is why we used powers of 2 here); it remains an open question whether the overhead of building 
that matrix outperforms (the additional) sorting. 
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A Appendix 


Proof of Lemma 2.1 Consider the lower bound for VaRQ,(L'’'). By De Morgan’s Law and Boole’s 
inequality, the distribution function F]^+ of L~^ satisfies 


d d d 

F:^,{x) = ¥{Y.L,<x) < F(min Lj < x/d) = [J {Lj < x/d}'^ < ^ IP(Rj < x/d) 
i=i ^ i=i i=i 


< dmaxFj{x/d). 


Now dmaxj Fj{x/d) < a if and only if x < dmiuj F~{a/d) and thus VaRc(-L''') = F“+(a) > 
d miuj F~ (a/d). 

Similarly, for the upper bound for VaRa(L'*“), we have that 

d 

Fl+{x) = F^ ^ Lj < > F(max Lj < x/d) = F(Li < x/d,L^ < x/d) 

S=i ^ 

d d 

= 1 - f(|J{Lj > x/d}) > max|l — o| 

i=i i=i 

d 

= max| ^ Fj{x/d) — d + 1, o| > maxjdmin Fj{x/d) — d + 1,0}. 

^ i=i ^ 


Now dmiuj Fj{x/d) — d+ 1 > a if and only if x > dmaxj F- ((d — 1 + a)/d) and thus VaRa(L^) = 
F~^{a) < dmaxj F~{{d — 1 + a)/d). □ 


Proof of Proposition \2.t^ 

1) Let s> s' and t' e [0, such that D{s',t') = D{s'). Define 

s — (s' — t'd) s — s' , 

*= a - —+ * 

so that 0 < t' < t < Let k = s' — t'd = s — td. If k > 0, noting that F is decreasing and that 
t' < t, we obtain 

d / _ \ 

D(s', t') — D{s, ~ J F{x)dx — J F{x)dxj > 0, 

so that D{s) < D{s,t) < D{s',t') = D{s'). If rc = 0 then D{s') = D{s', = dF{^) > dF{/i^) > 

D{s). 

2) Recall that D{s,t) = Fix) dx. Using the transformation z = (x — t)/is — td), we 

have 


D{s, t) = d [ F{sz + t{l — zd)) dz 

Jo 

Dehne C = {(s,f) | 0 < s < oo, 0 < t < and note that C is convex. Furthermore, if F is 
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convex, then D(s, t) is jointly convex in s and t on C since for A £ (0,1), 


D(y\s\ + (1 — A)s 2, Ati + (1 — A)t2) 

= d [ F{{Xsi + (1 - A)s 2 )- 2 : + (Ati + (1 - A)t 2 )(l - zd)) dz 
Jo 


< 


d [ F{X{siz + ti(l - zd)) + (1 - X){s2Z + t2{l 
Jo 

[ -^-^(((^1 + ^l(l ~ zd))) + (1 — A)-F((s2 + t2{l 
Jo 

XD{si,ti) + (1 - X)D{s2.,t2) 


— zd))) dz 

— zd)))) dz 


□ 


Proof of Proposition 2.6. First consider 0 7 ^ 1. Using (§, one can rewrite h{c) as 


h{c) = 


1 — a — dc 




Multiplying with c^^^d and rewriting the expression, one sees that h{c) = 0 is equivalent to hi{xc) = 0 
where 


Xc 


1 — a 
c 


(d-1) 


( 11 ) 


(which is in ( 1 , 00 ) for c £ ( 0 , (1 — a)/d)) and hi{x) = ^ ^x-i * -(('^ “ 1 )^ “ 1 “ !)■ H 

to see that hi{x) = 0 if and only if h 2 {x) = 0 , where 

h 2 {x) = {d/{l — 9) — — {d — l)x~^^^ + X — {d9/{l — 0) + 1), x £ (l,oo). (12) 


We are done for 0 7 ^ 1 if we show that /12 has a unique root on (l,oo). To this end, note that 
62 ( 1 ) = 0 and limj,.'|-oo h 2 {x) = 00 . Furthermore, 

/i2(x) = (1 — l/ 6 ){d/{l — 9) — + {d — l)9x~^^^~^ + 1, 

/i"(x) = {d + 9- - (1/0 + l)(rf - l)/9x-^/^-\ 

It is not difficult to check that h'f^x) = 0 if and only if x = (which is greater than 1 for 

d > 2). Hence, /12 can have at most one root. We are done if we find an xq £ (l,oo) such that 
h 2 {xo) < 0 , but this is guaranteed by the fact that limj,^! /i 2 (x) = 0 and limj,^! h 2 {x) = —{d— 2 ) /0 < 0 
for d > 2 . 

The proof for 9 = 1 works similarly; in this case, 62 is given by 

h 2 {x) = x^ + x(—dlog(x) + d — 2) — {d — 1), x £ (1, 00 ), (13) 


and the unique point of inflection of /12 is x = d/2. □ 

Proof of Proposition \2. Tj First consider q and 9^1. Instead of h, (0 and 0 allow us to study 

h 2 {x) = (d/(l — 9) — l)x“^/^’''^ — (d — l)x“^/® + X — (dd/(l — d) + 1), x £ [l,oo). 

Consider the two cases 9 £ (0,1) and 9 £ (1, 00 ) separately. If d £ (0,1), then d/(l —d) —1 > d—1 > 0 
and > x“^/^ for all x > 1, so /i 2 (x) > ((d/(l — d) — 1) — (d — l))x“^/^ + x — (dd/(l — d) + 1) > 

X — (dd/(l — d) + 1) which is 0 if and only if x = dd/(l — d) + 1. Setting this equal to Xc (defined 
in (111) and solving for c leads to the c; as provided. If d £ (l,oo), then using x“^/® < 1 leads to 
h 2 {x) > (d/(l — d) — l)x~^/^~^^ + X which is 0 for x > 1 if and only if x = (d/(d — 1) + 1)^. Setting 
this equal to Xc and solving for c leads to the q as provided. 
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Now consider 6 = 1. Similar as before, we can consider 0 and ( [I^ . By using that logx < 
and X > —for x £ [1, oo), we obtain h 2 {x) > x‘^+x{—dx^/^+d—2) — {d—l) > — 

which is 0 if and only if x = (d + Setting this equal to Xc and solving for c leads to the q 

as provided. 

Now consider Cu- It is easy to see that the infle ction point of /i 2 provides a lower bound Xc on the 

of inflection is x = Xc := 
for c then leads to c„ as stated. □ 


root of /i 2 . As derived in the proof of Proposition 2.0 tne point 
for 0 7 ^ 1 and x = d/2 for 6 = 1. Solving Xc = (1 — a)/c — {d— 1) 
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